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ABSTRACT 

Previous theoretical studies have found that repeating outbursts can occur 
in certain regions of an accretion disk, due to sudden transitions in time from 
gravitationally produced turbulence to magnetically produced turbulence. We 
analyze the disk evolution in a state diagram that plots the mass accretion rate 
versus disk surface density. We determine steady state accretion branches that 
involve gravitational and magnetic sources of turbulence. Using time- dependent 
numerical disk simulations, we show that cases having outbursts track along a 
nonsteady 'dead zone' branch and some steady state accretion branches. The 
outburst is the result of a rapid inter-branch transition. The gravo-magneto 
outbursts are then explained on this diagram as a limit cycle that is analogous 
to the well-known S-curve that has been applied to dwarf nova outbursts. The 
diagram and limit cycle provide a conceptual framework for understanding the 
nature of the outbursts that may occur in accretion disks of all scales, from 
circumplanetary to protoplanetary to AGN accretion disks. 

Subject headings: accretion, accretion disks — magnetohydrodynamics (MHD) 
- planets and satellites: formation — protoplanetary disks — stars: pre-main 
sequence — galaxies: nuclei 



1. Introduction 

Disk turbulence plays a key role in the outward transport of angular momentum that 
per mits disk accreti o n. Two main sourc es of disk turbulence are gravitational instabil- 
ity (iPaczynskil Il978t lLodato fc Ricd 120041) and m agnetic instability due to the magneto- 
rotational instability (MRI; iBalbus fc Hawley 1991 ). For the magnetic instability to operate, 
a critical level of ionization is required for the gas to couple strongly enough to the magnetic 
field. In certain situations, a nonturbulent, high density, cool, weakly ionized disk layer 
known as a 'dead zon e', can become gravitationally unstable, while remaining magnetically 
stable ( lGammidll996l ). The gravitational instability can lead to turbulent heating, resulting 



-2 - 



in increased ionization of the gas. The ionization can then trigger the magnetic instability, 
causing a much higher level of disk turbulence and accretion, an outburst. After an out- 
burst, the remaining disk gas cools and is repleni shed by accreting gas from larger radii (e.g. 
Armitage. Livio &: Pringldl200ll ; IZhu et al.l 120091 ) . The outburst then repeats at later times 
in an approximately periodic manner. Such behavior provides a possible model for FU Ori 
outbursts in young stars. 

The gravitational instability is thought to be a natural outcome of dead zones in a 
layered disk, regions where the disk i s nonturbulent near the disk mid-plane, but turbulent 
due to the MRI near the disk surface (jGammidll996l ). The surface disk turbulence occurs as 
a consequence of surface ionization by externa l sources of ionization, such as cosmic rays or 
X-rays, that then permit the MRI to operate (IGlassgold. Najita fc Igeall2004l ). In a layered 
region containing a dead zone, steady state accretion is generally not possible and the dead 
zone gains mass from some of the accretion flow near the disk surface. As mass builds up in 
the dead zone, it can become self-gravitating. As discussed above, the self-gravitating state 
is turbulent and can trigger a disk outburst. 

The purpose of this Letter is to describe the disk outburst cycle in terms of transitions 
between steady state configurations of an accretion disk. We determine the limit cycle for the 
outbursts in a diagram that plots the mass accretion rate M versus disk su rface density S. 
A similar approach was taken to explain the disk instability in dwarf novae ( IBath fc Pringle 
19821 ; iFaulkner. Lin &: Papaloizoulll983l ). The dwarf nova thermal- viscous instability can be 
understood by the S-shaped curve of steady state configurations in such a diagram. In this 
paper, we determine an analogous curve for outbursts in layered disks with gravitational and 
magnetic sources of turbulence. Our results do not conform to an S-curve because the middle 
portion of the S is missing. That is, unlike the S-curve case, there is range of accretion rates 
for which there are no steady state accretion configurations. 

In Section [2] we describe the equations for the disk model. In Section [3] we find steady 
state solutions to the disk equations and the resulting M versus £ curves. In Section 0] 
we describe the results of numerical simulations and show that the instability cycle can 
be understood in terms of transitions between the stable solutions in the S-M diagram. 
Section \5\ contains the discussion and conclusions. 



2. Disk Model 



The surface density e volution in a n accretion disk is determined by mass and angular 
momentum conservation (jPringld Il98ll ). The disk model we adopt is essentially the same 
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as Armitage et al. (2001). The main difference is that we explicitly include the effects of 
both MRI and self-gravity together in the disk evolution equations, while they treated these 
effects in separate equations. We have verified that the numerical results that they obtained 
can be recovered using the equations we describe below. 

We consider a disk in Keplerian rotation with angular velocity f2 = y^GM/R 3 around 
a mass M at radius R. The disk turbulence is modeled by the effective viscosities due to 
MRI and self-gravity. The surface density evolution equation is 



<9£ _ 3 d f i d 
~dt~RdR\ 2 dR 



'm^m ^g^g) 



(1) 



where v m is the kinematic viscosity in the magnetic layers and z/ g is the kinematic viscosity 
due to turbulence associated with self-gravity of the gas that may act outside the MRI layers. 
S m is the surface density in the MRI active layers and 

Sg = S — S m (2) 

is the surface density outside the MRI active layers. £ g can contribute to the mass flux if 
z/ g >0. 

We assume that the disk surface layers are ionized by external sources to a maximum 
disk density depth of S C rit/2 on the upper and lower disk surfaces. Therefore, the disk 
surface layers always contain MRI turbulence. At a given radius, we assume the disk is 
sufficiently thermally ionized for MRI to operate throughout the vertical extent of the disk 
if the mid-plane temperature is above some critical value T crit . Therefore, if either T c > T crit 
or if E < S cr it at some radius, then the disk is MRI turbulent (active) at all heights, where 
T c is the temperature at the disk mid-plane (central temperature). Otherwise, there is either 
a dead zone layer or the disk is self-gravitating in the presence of magnetic surface layers. 

The temperature at the disk mid-plane T c evolves according to the simplified energy 
equation 

dT c 2(Q+ - Q_) 



dt c p E 



(3) 



(jPringle. Verbunt fc Wadei Il986l ; ICannizzol Il993l ). The disk specific heat for temperatures 



around 10 K is c p = 2.77?.///, where 1Z is the gas constant and \i is the gas mean molecular 
weight. The local heating owing to viscous dissipation is taken as 

9 

Q+ = (^mS m + Z/g£ g ) . (4) 

o 

To determine the local cooling rate, we assume that each annulus of the disk radiates as a 
black body so that 

Q- = °T e \ (5) 
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where T e is the temperature at the surface of the disk and a is the Stefan-Boltzmann constant. 
The kinematic turbulent viscosity in the magnetic layer is taken to be 

r2 



n 



(6) 



where the sound speed is c m = ^/TZT m //j with temperature in the magnetic layer T m . The 
disk is self-gravitating if the Toomre parameter Q < Q CT a, where 



Q 



7T 



(7) 



and the sound speed at the disk mid-plane is given by c g = ^1ZT C / Jl, where we approximate 
the temperature of the self-gravitating region that extends to the disk mid-plane as ~ T c . The 
effective kinematic viscosity from the turbulence induced by the self-gravitational instability 
is approximated by 



Qc 



Q 
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for Q < Qcrit and zero otherwise (ILin &: Pringldll987l . Il99dl ). 

The mid-plane disk temperature for an optically thick disk in thermal equilibrium is 
obtained by considering the energy balance in a layered model above the disk mid-plane. 
One layer contains the surface density S m /2. The other layer contains the complementary 
surface density £ g /2. The results are that 



4 9 2 
&T C = — f2 (V m S m T m + z/gEgT) 



and 



rpA 

m 



T T 



The optical depth to the magnetic region is 



3 . . S m 

and the optical depth within the complementary region is 



(9) 
(10) 

(11) 



with 



'12) 



(13) 
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Note r g is defined even in a dead zone layer (y z = 0) with E g defined by equation (j2J). We 
adopt the simplified opacity of Armitage et al. (2001) 

k(T) = 0.02 T°- s cm 2 /g. (14) 

The energy equation (j3J) in a steady state has the solution 

9 

aT e = -Q 2 (z/ m £ m + z/ g S g ) . (15) 

o 

From equations (l9l- (fT5l) . we obtain an expression for the cooling function 



Q_ = aT c 4 = r- 1 (aT c 4 + V^£ m T g J . 



(16) 



We apply this cooling function to equation ([3]), even when the disk is not in thermal equi- 
librium. This means that we do not attempt to treat the cooling during the transitions 
consistently. 



3. Steady State Disks and State Transitions 

We consider the disk to be supplied by material at a constant mass input rate into an 
outer region. In a steady state, the mass flux through the disk is constant in radius interior 
to the mass input region. We determine the relationship between the mass flux M and the 
disk surface density S at some radius. 

In a steady state, far from radius of the central object, the surface density equation ([I]) 
has the solution 

M = 3tt (z/ m £ m + z/gSg) . (17) 
From equations (fT5l) and ffTTl) . we find the steady state surface temperature of the disk 

T e = -= • 18 

\ QUO J 

Notice that T e is a function of the total accretion rate M. 

We sketch typical solutions, in the S-M plane at a given radius in Fig. [1] to illustrate 
their principal properties. We show the various branches of disk solutions, as described 
below. In Section HI we consider particular numerical solutions. 
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3.1. MRI Disk Branches 

If T c > T crit or S < S crit at a given radius, then the disk is assumed to be fully MRI 
turbulent, that is magnetically turbulent at all heights (E g = 0). The disk surface density 
in this case simply follows from equation f fT7|) 

M = 3vrz/ m S m . (19) 

At a fixed radius, equations flBJ), (fTUj) . f TTgj) . and (fT9~j) imply that M oc £m . 

3.2. Gravo-Magneto(GM) Disk Branch 

Over a certain range of accretion rates, there exists a branch of steady state solutions 
where the disk is gravitationally unstable and has MRI active surface layers with S m = £ CT it- 
We numerically determine M(S) over a range of £ for which Q < Q cvl t and T < T crit . 
Equations ^ and f|T7|) are numerically solved together with auxiluary equations ([HD - ©, 
(HDD - (HID, (USD, and ([IBD to obtain this function. 

3.3. Dead Zone Branch 

Between the lower MRI and the GM steady state branches there are disk configurations 
that involve dead zones. The dead zone is a nonturbulent mid-plane layer that lies vertically 
between MRI actively accreting surface layers. Unlike the other branches we have considered, 
this branch is not in steady state and therefore not described by steady state solutions to the 
disk equations. It is plotted as a dashed line in lower portion of Fig. [TJ If a disk lies along 
this line, it will evolve to the right, as the dead zone gains mass, while the disk accretion 
rate is pinned to the rate provided by MRI turbulence that resides in surface layers of fixed 
density E crit . The mass gain of the dead zone is produced by some of the accretion flow in 
the disk surface layers that becomes incorporated into the dead zone. Over time the dead 
zone gains sufficient mass for the disk to become self-gravitating. The disk state then enters 
the GM branch. 



3.4. Accretion rates with no steady state 

There exists a range of mass accretion rates for which no steady state solution exists. 
This range is indicated by the lower shaded region in Fig. [TJ This gap comes about because 
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Fig. 1. — A schematic diagram of steady disk solutions in the S-M plane at some radius in 
the disk. The solid lines show steady state solutions. The two lower straight line branches 
on this log-log plot show the disk configurations with magnetic (MRI) turbulence at all 
heights. The lower branch ends where E = £ m = S cr j t . The higher branch begins where 
T c = T cri t- The curved GM branch is for a disk configuration with magnetically driven surface 
turbulence that overlies a mid-plane region that is turbulent due to the effects self-gravity. 
The dashed line corresponds to nonsteady configurations in which a dead zone (nonturbulent 
region) underlies surface layers with MRI turbulence. The lower shaded region contains no 
steady state solutions for any value of M. The upper S-curve describes the thermal-viscous 
instability. The dotted line shows its unstable steady state solutions. 
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the GM branch terminates at an M for which T c = T cr ; t , while the higher MRI branch starts 
at a larger value of M that is also at T c = T crit . As we will see in Section HI a disk that has 
an accretion rate lying in this range will quickly transition to either of these steady state 
branches. As a consequence of such transitions, the disk undergoes outbursts. 

3.5. Global Steady State 

We are considering a disk that undergoes accretion from an external source at a constant 
accrete rate. The E-M diagram applies at each disk radius. The range of unstable accretion 
rates shown as the lower shaded region in Fig. [T] varies with radius. For a disk to be globally 
stable against gravo-magneto outbursts, there should be no radius at which the disk accretion 
rate lies within this range. 

3.6. Thermal-Viscous Instability 

The case of outbursts involving the thermal-viscous instability occurs at higher disk 
temperatures and accretion rates than the gravo-magneto instability. This regime is sketched 
as the well-known S-curve in the upper portion of Fig. HJ The instability occurs along the 
middle portion of the S-curve (the upper dotted line within the upper shaded region). This 
middle portion consists of unstable steady state solutions that occur between other the two 
stable branches. The situation with the gravo-magneto outbursts is different because there 
are no intermediate unstable steady state solutions, i.e., there are no steady state solutions 
in the lower shaded region. For a given accretion rate, it is possible for the disk to develop 
the thermal- viscous instability at small radii and the gravo-magneto instability further out. 

4. Time Evolution 

We analyze the time evolution of a disk subject to only the gravo-magnetic instability. 
We consider a case with R = 3 AU, M = 1M , M = lO~ 6 M yr _1 , Q crit = 2, E crit = 
200g/cm 2 , T cri t = 800 K, /i = 2.3, and a m = a g = 0.01. The steady state solutions plotted 
in Fig. [2] show that no steady state branch exists at this mass accretion rate. Therefore, we 
expect outbursts in the accretion rate on to the central star. 



We time- integrate the evolution equations (pQ) and (j3J), together with equations and 
([IE]) , for a dis k that extends from R\ n = 2.33 x 10~ 3 AU to R nut = 40 AU around a solar mass 
star (see e.g. lArmitage. Livio fc Pringldl200ll ; iMartin fc Lubowll201ll ). In this calculation, 
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lo gio( E /g cm- 2 ) 

Fig. 2. — The E-M s and T e plane for a disk parameters listed in SectionlU M s = 37r(z/ m S m + 
z/ g Sg) is the steady state accretion rate that applies to the steady state branches. The dashed 
line and thick lines disk solutions are defined in Fig. [TJ The thin line shows the limit cycle 
of time-dependent disk evolution for T e versus S for an accretion rate on to the disk of 
M = 10 -6 M Q yr _1 for which no steady state branch exists having M = M s . 
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we ignore the effects of the thermal-viscous instability shown by the S curve in the upper 
part of Fig. [TJ The grid consists of 120 points distributed uniformly in R^. Material is 
added to the disk at a constant rate of M — 10 -6 M yr _1 at a radius of 35 AU. At the inner 
radius, we impose a zero torque boundary condition, so that there is an inward flow of gas 
out of the grid and towards the central star. The flow is prevented from leaving the outer 
boundary by means of a zero radial velocity outer boundary condition. 

The disk undergoes outbursts that reach a limit cycle after a few outbursts. In Fig. [2] 
we show the disk evolution over a limit cycle. The cycle proceeds in a counter-clockwise 
sense on the diagram. Along the low horizontal section of the cycle, the disk is in the dead 
zone state, as discussed in Section 13.31 The surface density increases and the disk evolves 
to the right. The disk then becomes self-gravitating with Q < Q crit as a consequence of this 
accumulated dead zone mass. The disk state then moves along the gravo-magneto branch 
until it reaches the critical temperature for the MRI to act, T cr i t . The disk then becomes 
well enough ionized for MRI to act and the disk state jumps up to the fully active MRI 
branch. The upward jump corresponds to the start of a disk outburst that proceeds at a 
higher accretion rate. The mass then drains out of the disk and the disk cools. As a result, 
the path then moves down and to the left along the MRI branch. Once the temperature T c 
falls below the critical value T cr it, the path moves downward. The disk once again forms a 
dead zone and the process repeats itself. 

5. Discussion and Conclusions 

We have explained the gravo-magneto outbursts triggered in a layered accretion disk by 
means of a E-M diagram, along the lines of the S-curve explanation for the thermal-viscous 
instability in dwarf nova outbursts. We determined steady state disk solutions that compose 
branches in the diagram. At a given radius, there are two types of steady-state solutions (see 
Fig. [1]). One type, labelled MRI, involves a disk that is magnetically unstable and turbulent 
at all heights. This type describes two branches of solutions, one at low and the other at 
high accretion rates. Another type, labelled GM, is gravitationally unstable and turbulent 
near the mid-plane and is magnetically unstable and turbulent only near the disk surface. 
This type describes one branch of solutions. At a given radius, there may be a range of 
disk accretion rates for which there are no steady state solutions, the lower shaded region in 
Fig. [TJ This range of accretion rates varies with radius. If there is a radial zone within a disk 
that has an accretion rate lying in this range, then the disk will undergo outbursts in that 
zone. The outbursts are understood in terms of a limit cycle in the M(E) diagram as shown 
in Fig. [2j The cycle tracks along steady state branches of the disk solutions and along the 
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'dead zone' branch. The cycle jumps between states of low and high accretion rates at the 
initiation and termination of an outburst. The outburst model may be applied to accretion 
disks of all size scales from circumplanetary to protoplanetary to AGN disks. 

We have made several simplifications in our model that will be investigated further in 
future work. We have ignored the thermal-viscous instability that produces the dwarf nova 
outbursts that is sketched in the upper portion of Fig. [TJ We adopted a simplified opacity 
that does not produce this curve in the numerical models. The thermal-viscous instability 
could be triggered at smaller radii than the GM instability. The interaction between the two 
instabilities may produce some interesting behavior applicable to T Ta uri stars. The role 
of t he thermal - yiscou s instability and its S-curve has been considered by iBell fc Lin! (119941 ) 
and Izhu et all (I2009I ). 



Following previous work (Armitage et al. 2001), we have taken a rather small value for 
the viscosity in the active regions of a m = 0.01 and a large value for the depth of the surface 
magnetic turbulence of £ r .rit = 200g/cm 2 . Some general considerations suggest that a ~ 0.1 
(IKing. Pringle fc Livid 120071 ) . We note that if the a m were higher, then the disk would be 
more stable against the gravo-magneto outbursts. Similar results to those plotted in Fig. [2] 
are obtained for a m = 0.1, but with E crit = 20g/cm 2 . 



We have assumed that the turbulent viscosity associated with self-gravity, a g , depends 
on the Toomre Q parameter. However, it has been shown in several papers that a 3 actually 
depends on the cooling rate (jGammidl200ll ; ICossins. Lodato &: Clarkd 120091 ). The final disk 
Q value is found to be self-regulated at a value ~ 2. 

We have neglected external heating of the disk and assumed that the accretional heating 
dominates. For the parameters associated with the GM instability in Fig 2, this assumption 
may be valid. However, more generally external heating should be considered. 

The use of a fixed value for the maximum depth of nonthermal ionization of the 
outer magnetic layers, Ecrit) is an approximation. A more accurate approach to apply 
a critical magnetic Reynolds number required for the MRI to operate, along the lines of 
Matsumura &: Pudrita (120031 ). In addition, we have taken a single critical temperature for 
thermal ionization T cr i t at the disk mid-plane temperature above which the disk becomes 
MRI turbulent at all heights. Instead, there may be a thin turbulent layer near the disk 
mid-plane tha t develops at this mid-plane temp erature that increases in thickness at higher 
temperatures (IZhu. Hartmann fc Gammidl2009i ). 



Several authors have suggested that there is a small but non-zero turbulent viscosity 
in the dead zone that develop s as a response to tur bulence driven in the magnetic surface 
layers ([Fleming fc Stone! 120031 ; iTurner fc Sand 120081 ). The range of steady state disk flow 
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solutions is then increased (e.g., iTerqueml 120081 ) . To some extent, this effect will stabilize 



the outbursts as material can flow through the (nearly) dead zone. However, there may be 
issues with the model if the rate of flow through the dead zone approaches that through the 
magnetic surface layer. Such topics can be explored within the framework we have described. 
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